Three-dimensional topological phases in a layered honeycomb spin-orbital model 
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We present an exactly solvable spin-orbital model based on the Gamma-matrix generalization of 
a Kitaev-type Hamiltonian. In the presence of small magnetic fields, the model exhibits a critical 
phase with a spectrum characterized by topologically protected Fermi points. Upon increasing 
the magnetic field, Fermi points carrying opposite topological charges move toward each other and 
annihilate at a critical field, signaling a phase transition into a gapped phase with trivial topology 
in three dimensions. On the other hand, by subjecting the system to a staggered magnetic field, an 
effective time-reversal symmetry essential to the existence of three-dimensional topological insulators 
is restored in the auxiliary free fermion problem. The nontrivial topology of the gapped ground state 
is characterized by an integer winding number and manifests itself through the appearance of gapless 
Majorana fermions confined to the two-dimensional surface of a finite system. 



I. INTRODUCTION 

Topological phases of matter are one of the most 
remarkable discoveries in modern condensed-matter 
physics. Instead of a local order parameter describ- 
ing broken symmetries of the system, this novel state of 
matter is characterized by a topological invariant which is 
insensitive to small adiabatic deformations of the Hamil- 
tonian. A classic example is the integer quantum Hall 
(IQH) effect, where the quantized Hall conductance cor- 
responds to a topological invariant called the first Chern 
number or the TKNN integer.^ Another distinctive prop- 
erty of IQH states is the appearance of zero-energy modes 
at the sample edge, despite the fact that all bulk excita- 
tions are fully gapped. Due to their topological nature, 
these conducting edge modes persist even in the pres- 
ence of disorders. Recently, theoretical investigations 
have shown that similar topological insulators can be 
generalized to time-reversal invariant systems (quantum 
spin-Hall effect)®"^ and to three dimensions.^ -"^^ Signa- 
tures of topological insulators such as quantized con- 
ductance and protected surface Dirac cone have been 
reported experimentally in semiconducting alloys and 
quantum wellsj^^— Following these developments, sys- 
tematic classifications of topological insulators have also 
been proposed.—"— 

Although most discussions of topological insulators are 
in the context of tight-binding fermionic models or mean- 
field superconductors, it has been shown that topological 
insulators can also emerge from strongly correlated elec- 
tronic systems i^"— An exactly solvable example is Ki- 
taev's anisotropic spin-1/2 model on the two-dimensional 
honeycomb lattice.— As shown in his seminal paper, the 
spin model can be reduced to a problem of free Majo- 
rana fermions coupled to a static Z2 gauge field. The 
ground state of the Kitaev model has two distinct phases. 
The gapped Abelian phase is equivalent to the toric code 
model^- whose excitations are Abelian anyons. Relevant 
to our discussion is the non- Abelian B phase in the pres- 
ence of a magnetic field. This phase is characterized by 
an integer winding number v = ±1. Similar to IQH 
insulators, the non-Abelian phase also supports gapless 



chiral edge modes (whose chirality depends on the sign of 
magnetic fields) except that the edge modes here are real 
Majorana fermions, as contrasted to complex fermions in 
the case of IQH states. 

There has been much effort to generalize Kitaev 
model to other trivalent lattice a^^'^^ and to three- 
dimensions j^"— Probably the most notable example is 
the discovery of a chiral spin liquid as the ground state 
of Kitaev model on a decorated-honeycomb lattice;^ On 
the other hand, despite being exactly solvable, we find 
that most gapped phases of 3D Kitaev model is topo- 
logically trivial.-- The fact that the model can only be 
defined on lattices with coordination number 3 signifi- 
cantly constrains the possible Majorana hopping Hamil- 
tonian. Recently, noting that the exact solvability of 
Kitaev model relies on the fact that the three spin-1/2 
Pauli matrices realize the simplest (dimension-2) Clifford 
algebra, the so-called P-matrix generalization of the Ki- 
taev model offers richer possibilities of engineering ex- 
actly solvable models with unusual phases in both 2D and 
3D.—"— Physically, models based on, e.g. dimension-4 P 
matrices can be interpreted as spin-| models, spin-^i 
models, or spin-orbital models. 

Based on the above P-matrix generalization of Kitaev 
model, emergent topological insulators have been demon- 
strated on a 3D diamond lattice with coordination num- 
ber iM-i^ To construct a Kitaev-type model on such a 
lattice, one needs four mutually anticommuting opera- 
tors for the four nearest-neighbor links. This can be 
easily realized using the dimension-4 P-matrix represen- 
tation of the Clifford algebra. By exploiting the redun- 
dancy of representing two such sets of bosonic operators 
in terms of 6 Majorana fermions, the diamond-lattice 
model enjoys an effective time-reversal symmetry (TRS) 
which is essential to the existence of topological insu- 
lators in 3Dii^ With only nearest-neighbor interactions, 
the diamond-lattice model reduces to a problem of two 
identical copies of free Majorana fermions sharing the 
same Z2 gauge field. The permutation symmetry be- 
tween the two fermion species thus manifests as a TRS. 
In the weak-pairing regime of the model, hybridization 
of the two Majorana species results in a gapped ground 
state with nontrivial topology. 



A natural question to ask is whether it is possible to re- 
alize the required TRS in a Kitaev-type model via generic 
physical mechanisms. In this paper, we provide such an 
example through a natural generalization of the honey- 
comb Kitaev model. Using a similar F-matrix formalism, 
we study a Kugel-Khomskii-type^^ spin-orbital Hamilto- 
nian on a layered honeycomb lattice whose coordination 
number is 5. We further consider perturbations due to 
a weak magnetic field and single-ion spin-orbit interac- 
tion. In the weak-pairing regime of our model, the ground 
state remains gapless up to a critical field strength. The 
Fermi points of this critical phase are characterized by 
a nonzero topological invariant, hence are stable against 
weak perturbations. As the field strength is increased, 
Fermi points with opposite winding numbers eventually 
annihilate with each other and the spectrum acquires a 
gap above a critical field. Since the TRS is explicitly 
broken in the presence of a uniform magnetic field, the 
gapped phase represents a multilayer generalization of 
the 2D Kitaev model, and is topologically trivial in three 
dimensions. 

On the other hand, when the sign of magnetic field 
alternates between successive honeycomb planes, an ef- 
fective TRS is restored in the auxiliary fermionic model, 
and the corresponding ground state in the weak-pairing 
phase is equivalent to a topological superconductor be- 
longing to symmetry class Dili in Altland-Zirnbauer's 
classification4S We also show that the quantum ground 
state is characterized by an integer winding number, con- 
sistent with the general classification of 3D topological 
insulators A physical consequence of the nontrivial 
ground-state topology is the appearance of surface Ma- 
jorana fermions which remain gapless against perturba- 
tions respecting the discrete symmetries of the Hamil- 
tonian. The specific multilayer geometry of our model 
also allows us to analytically demonstrate the existence 
of surface Majorana fermions obeying a Dirac-like Hamil- 
tonian by relating the surface modes to the chiral edge 
modes of 2D Kitaev model. 



II. THE SPIN-ORBITAL MODEL 

We study a Kugel-Khomskii-type^ spin-orbital model 
defined on a layered honeycomb lattice shown in Fig. [T] 
In addition to a spin- 1/2 degree of freedom, each lattice 
site has an extra doublet orbital degeneracy. The spin 
and orbital (pseudospin) operators are denoted by Pauli 
matrices cr" and t" (a — x,y,z), respectively. As in the 
original honeycomb Kitaev model, nearest-neighbor links 
lying on a honeycomb plane are divided into three types: 
X, y, and z, depending on their orientations. The model 
Hamiltonian is defined as follows: 
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FIG. 1: A layered honeycomb lattice. The five distinct 
nearest-neighbor links are indicated by the corresponding ex- 
change constant. The open and filled circles denote sites be- 
longing to the two inequivalent sites of a unit cell. The primi- 
tive vectors of the underlying Bravais lattice are ai — (1, 0, 0), 

a2 = (|, -^,0), and as = (0,0, 1). An arrow from site j to 
site k means the corresponding link variable Ujk ~ +1- 



Here j runs over the lattice sites, j ± z denote nearest 
neighbors along the two vertical links, and j + Sa de- 
notes in-plane nearest neighbor along the link of type a. 
The prime in the first term indicates that the summa- 
tion runs over sites on every second honeycomb layer. 
The exchange constant is Jy on vertical links, and is Ja 
on a- links lying on a honeycomb plane (see Fig.[T]). The 
spin-orbital interaction within each honeycomb layer re- 
sembles the original 2D Kitaev model with spin-1/2 oper- 
ator cr" replaced by the spin-orbital operator t^ct". The 
inter-layer interaction in our model, on the other hand, 
involves only orbital operators; the orbital interaction al- 
ternates between the and types along successive 
vertical links. In addition to discrete lattice symmetries, 
the Hamiltonian is invariant under a tt rotation about 
axis followed by lattice translations along the z direction 
by one layer. 

Since the five spin-orbital operators r^, t^, and t^ct" 
anticommute with each other, they generate a 4 x 4 
matrix representation of the Clifford algebra. With an 
enlarged Hilbert space, one may introduce 6 Majorana 
fermions c and 6^ (/i = 1, • • • , 5) such that: 

= ib^c, T^'oy = ib^c, T^'a'' = ib^c, 

= ih^'c, ry = ib^c. (2) 

By denoting these operators as F^ = ifo^c, Hamilto- 
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FIG. 2: Configuration of Z2 gauge fields {ujfc} and {mj^} on 
(a) square and (b) hexagonal plaquettes; Ujk = 1 if there is 
an arrow pointing from site j to site k. The numbers 1 and 2 
label the two inequivalent sites within a unit cell. 

nian ([!]) can be recast into a Kitaev-typc interaction 

/I— 1 /I— links 

with Ji ~ Jxi Ji — Jyi a-nd J3 = Jz for links lying on a 
honeycomb plane, and J4 = J5 = Jy along the vertical 
links. The six Majorana fermions form an 8-dimensional 
Hilbert space which is twice as large as the local physi- 
cal Hilbert space. This redundancy can be remedied by 
demanding the allowed physical states be eigenstate of 
gauge operator D = ic n^=i with eigenvalue +1. This 
is also consistent with the identity r^r'^r^a'^aya^ = — 1. 
We note that the same F-matrix extension of the Ki- 
taev model on a 2D decorated square lattice (Shastry- 
Sutherland lattice) has been studied in Ref. [H. 

Following Kitaev)^ we introduce the link operator 
Ujk = ibjb'j^, where fi = fj,jk implicitly depends on the 
type of link connecting sites j and k. The Hamiltonian 
then becomes 

■Ho = J^UjfcCjCfc. (4) 

jk 

A remarkable feature of the fermionic Hamiltonian first 
noted by Kiatev is that the link operators ujk commute 
with each other and with the Hamiltonian. We may thus 
replace them by their eigenvalues ujk = ±1, which act 
as an emergent Z2 gauge field. Consequently, for a given 
choice of {ujfc}, Hamiltonian Q reduces to a problem of 
free Majorana fermions with nearest-neighbor hopping 
ijk = JfiUjk- However, since Ujk does not commute with 
operator Dj , the spectrum of the free fermion Hamilto- 
nian depends only on gauge invariant quantities which 
are given by the product of link operators around the 
boundary of elementary plaquette: Wp = Y[{jk)£dp''^3''- 
As Wp = 1, the flux associated with a given plaquette is 
also given by a Z2 variable Wp = ±1. 

The question remains of what choices of Ujk, or equiv- 
alently the Z2 fluxes Wp, give the lowest ground-state 
energy of the free fermion problem. Numerically, we find 
that such a state is attained when all hexagons are vortex 
free while all squares contain a tt flux: 

WQ = 1, wa = -1, (5) 



consistent with Lieb's theorem^S A gauge choice which 
gives the desired plaquette fluxes without enlarging the 
unit cell is shown in Fig. [21 With this specific choice 
of Z2 gauge field, the fermion spectrum can be obtained 
analytically using Fourier transformation. 

In the weak-pairing regime to be discussed below, 
model (dl) exhibits a gapless phase similar to the B phase 
of 2D Kitaev model. In order to explore possible topo- 
logical insulators emerging from this critical phase, we 
consider the following single-ion perturbations: 

j J 

The first term is the Zeeman coupling to an external 
magnetic field hj = Ty^h, where the phase factor 77^- = ±1 
depends only on the z coordinate of spin j , i.e. the field is 
constant within a given honeycomb plane. As mentioned 
previously, we shall consider two special cases in this pa- 
per: uniform field with rjj — 1 and staggered field with 
rjj — (— . Since spins transform as Ta'^T~^ — — cr" 
under time-reversal T, the first term above explicitly 
breaks TRS. 

The second term in Eq. ^ represents a spin-orbit-like 
interaction where A is an effective coupling constant, and 
n is a unit vector specifying the local anisotropy axis. 
Such a coupling arises when the orbital basis = ±1) 
carries a nonzero orbital angular momentum L ~ r^fi 
which is parallel or antiparallel to the local anisotropy 
axis depending on = +1 or —1, respectively. An ex- 
plicit example is given by two degenerate orbitals with 
t2g symmetry, e.g. \yz) and \zx) in a tetragonal crystal 
field. By identiiying |t^ = ±1) as \yz) zLi\zx), respec- 
tively, the orbital basis has a nonzero angular momentum 
pointing along the symmetry axis of the tetragonal field. 

The effects of Hi can be studied following the pertur- 
bation treatments discussed in Ref. [2^. Essentially, one 
constructs an effective Hamiltonian acting on the sub- 
space which is free of vortex-type excitations. The first- 
order correction vanishes identically as both single-ion 
perturbations create tt fiuxes on hexagons^^ whereas the 
second-order terms simply modify the nearest-neighbor 
exchange constants Ja- Nontrivial corrections to the 
fermion spectrum arise at the third-order perturbation 
which involves multiple spin-orbital interactions, e.g. 

Such terms introduce an effective second nearest- 
neighbor hopping for Majorana fermions 

jk 

where k ^ X^h/J'^ and 77 = ±1 is a constant within the 
plane containing sites j and k. The additional second- 
neighbor Z2 field u'jf, — —u'f,j is shown by the dashed 
line in Fig. [2j Depending on the direction of n and h, 
the hopping amplitude k could take different values along 
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inequivalent second-neighbor links. As the main purpose 
of this term is to introduce a spectral gap, to avoid un- 
necessary complications, we shall assume the symmetric 
case in the following discussion. 



III. UNIFORM MAGNETIC FIELD 

We first discuss the ground state in the presence of a 
uniform magnetic field, i.e. rjj = +1 for all layers. As 
will be discussed below, the gapped phase at large fields 
is characterized by a tt2 Chern number corresponding to 
a multi-layer generalization of the 2D Kitaev model, and 
is topologically trivial in 3D (which generally is charac- 
terized by TTa homotopy groups). This is mainly because 
the TRS essential to the existence of 3D topological insu- 
lators is explicitly broken by the uniform field. However, 
the critical phase in the case of small fields is interest- 
ing in itself and is similar to the gapless A-phase of '^He 
discussed in Ref. [1; the gaplessness of both phases are 
topologically protected. 

Since the original unit cell of the lattice is preserved 
by the Z2 fields ujk and u'jf. shown in Fig. [2J we ex- 
press the site index as j = (r, s), where r denotes the 
position of the unit cell, and s = 1,2 indicates the two 
inequivalent sites in a unit cell. With Fourier transforma- 
tion ak,s = X^r e"*'^'^'""'"'*^^ /V'2N, where dg is a basis 
vector and N is the number of unit cells, the Hamilto- 
nian becomes % = + Hi = \ Y.\, *k-f^(k) *k, with 
*k = (ak,i, ak,2)^ and 



g(fc,)-|-A(k^) 
»/(k±)* 



.g(fc,)-A(k^) )' 



i/(k) ^ 

= Im/(ki)r"+Rc/(ki)T^+ [g(fc,)-fA(ki)]r 
For convenience, we have defined the following functions: 



-l{k,) = 4J|| sinfc^, /(k_L) = 2 ^ J„ e' 

o 

A(kj^) = 8k sin ^ (^os 



k I (So 



(9) 
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The three vectors ^1^2 = (^'5^^)' ^'^'^ ^3 = (0' 75) 
connect nearest-neighbor sites in a honeycomb layer. In 
the following we shall focus on the emergent free fermion 
problem. The Pauli matrices t'^ appearing in the single- 
particle Hamiltonian now act on the sublattice index s 
(not to be confused with orbital pseudospins). 

The hermitian matrix (k) in Eq. also satisfies 



H^{-k) = -H{k), 



(10) 



which is the defining property of symmetry class D in 
Altland-Zirnbauer's classification;^ As 3D insulators in 
this class is topologically trivial)^ the gapped phase of 
Eq. ([5]) represents a trivial multilayer generalization of 
the 2D Kitaev model. Nonetheless, for small magnetic 



fields such that the second-neighbor hopping k ^ X^h is 
below a critical value Kc, the fermion spectrum remains 
gapless in the weak-paring regime of the model. The cor- 
responding critical phase is characterized by topologically 
protected Fermi points as we shall discuss below. 
Diagonalizing Hamiltonian ([5]) yields a spectrum: 



e±(k) = ±V [g{h) + A(k^)] + |/(k^)p. (11) 

When one of the in-plane coupling is much larger than 
the other two, e.g. Jz ^ Jx,Jy, there is no solu- 
tion to equation f{k.±) = 0, and the spectrum is al- 
ways gapped irrespective of the applied magnetic field. 
This phase corresponds to the A phase of the 2D Kitaev 
model, ^® and is referred to in the following as the strong- 
pairing phase based on analogy with the p-wave topolog- 
ical superconductors 4i On the other hand, in the weak- 
pairing regime of the model where the three in-plane cou- 
plings satisfy the triangle inequalitieSj^^ the model dis- 
plays a possible gapless phase as two solutions k^ exist 
for the equation f(k±) = 0. To be specific, we now con- 
centrate on the symmetric case Jx = Jy — Jz = J, where 
zeros of /(kj_) are at the corners of the 2D hexagonal 
BriUouin zone k*j_ = {±^,0). 

In contrast to the 2D Kitaev model where applying a 
magnetic field immediately opens an energy gapj^ spec- 
trum (|11|) of the 3D model remains gapless when k is less 
than a critical strength \k\ < Kc. For symmetric in-plane 
couplings, we find Kc = 2J||/3a/3 and the nodes of the 
spectrum e±(k) are located at 
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0. 



(12) 



where ^ — arcsin(K/Kc). These Fermi points are topo- 
logically protected and are robust against weak pertur- 
bations. To see this, we first rewrite Hamiltonian ([5]) 
as H{k) — e_(-(k) rh(k) -r, where m(k) is a unit vector. 
A topological invariant characterizing the singularities of 
the spectrum is given by the winding number of map- 
pings from a sphere enclosing the Fermi point k* to 
the 2-sphere of the unit vector rhi^ 



^e^'^P m • (9^m x dura) 



(13) 



The above definition corresponds to the second homotopy 
group TT2{S'^) = Z, which characterizes point defects in 
an 0(3) spin field4^ 

Examples of topologically nontrivial Fermi points are 
given by the spectra of Weyl Hamiltonian describing a 
massless spin-1/2 particle: -ffwoyi = HcT^d^, where c 
is the speed of light. ^ The plus and minus signs refer 
to left and right-handed particles, respectively. Using 
Eq. (jl3p the winding number of Weyl spinors can be read- 
ily computed, resulting v — ±1 for right and left-handed 
particles, respectively. Take left-handed Weyl spinor for 
example, the expectation value of its spin is parallel to 
its momentum: (t) || p. In the ground state with filled 
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FIG. 3: (Color online) (a) Projection of ground-state pseu- 
dospin (r) — — m(k) on a face of the Brillouin zone. The two 
vertical edges of the face correspond to K = {^,0,kz) and 
K' = {^,^^,k^) = (-^,0,fcz). The location of the four 
topologically non-trivial Fermi points are indicated by black 
arrows, (b) Dispersion along the K' edge of Brillouin zone 
for various k ~ X^h. Upon increasing magnetic field h, the 
two Fermi points carrying opposite winding numbers move 
in opposite directions along the K' edge. At a critical field 
(k = Kc), the two singularities merge to form a new Fermi 
point with trivial winding number v = at { — ^, 0, ^). The 
spectrum becomes gapped above Kc- 



negative-energy states, Fermi point with 1^ = 1 thus looks 
like a magnetic monopole (a hedgehog) in momentum 
space. These singular points are robust in the sense that 
it is impossible to continuously deform a hedgehog into a 
uniform spin configuration (corresponding to the trivial 
case of v — 0). 

To compute the winding number of Fermi points in 
our 3D model, we expand Hamiltonian (|8]) around the 
singular points, e.g. 

iJ(k^ + p) = -V^PyT'' - V^p^T^ T (14) 

where the Fermi velocities are vj_ — and vn — 



The dispersion around these points 

^±v/KpIFTT^^- Af- 



has a conic singularity: e± '■ 
ter rotation and mirror-inversion, Eq. (jl4p is essentially 
equivalent to the momentum-space Weyl Hamiltonian 
discussed above. The four Fermi points form two 
pairs (labeled by r = 1,2) of singularities with opposite 
winding numbers — ±1. The ground-state configura- 
tion of pseudospin (r) = — m(k) projected onto a face of 
the Brillouin zone is shown in Fig.lSja). The two pairs of 
singularities are located at the two inequivalent edges K 
and K' of the 3D Brillouin zone. Instead of a monopole- 
like configuration, vector field m(k) in the vicinity of 
these Fermi points has a saddle-point-like singularity. 

Although these Fermi points are topologically pro- 
tected due to their nonzero winding numbers, a spectral 
gap can still be induced through mutual annihilation of 
Fermi points carrying opposite topological charges. This 
process is illustrated in Fig. [3] Take for example the two 
Fermi points on the iC'-edge of the Brillouin zone. Upon 
increasing magnetic field /i, the two singularities charac- 
terized by winding numbers v — ±1 move toward each 
other along the K' edge, and eventually merge to form 
a new Fermi point when k reaches the critical Kc- Be- 
cause of the conservation of topological charge, the new 
Fermi point has a winding number v — 0, hence is topo- 
logically trivial. Above the critical field, the spectrum is 
fully gapped. 

It is also interesting to understand the critical B phase 
of Kitaev's original 2D model from the perspective of 
topological winding number. As discussed in Ref. [2^, 
the gaplessness of B phase is protected by TRS on the 
bipartite honeycomb lattice. This discrete symmetry es- 
sentially forces the unit vector rh(kj^) to lie in the xy 
plane as terms in Eq. ^ is not allowed by TRS. For 
a planar unit vector m whose tip lies on a circle S^, the 
singularities are characterized by an integer topological 
invariant, also known as the vortex winding number, cor- 
responding to 7ri(S'^) = Zi^ The two Fermi poins in the 
B phase of Kitaev's model can be viewed as vortices car- 
rying opposite winding numbers = ±1, respectively. In 
the presence of perturbations breaking the TRS, the vec- 
tor rh now lives on a 2-sphere. Since 7Ti{S'^) — 0, singu- 
larities of the spectrum are then topologically trivial (7r2 
characterization of 2D singularities is not well defined). 

The gapped phase above the critical Kc represents a 
trivial generalization of the 2D Kitaev model in much 
the same way as the multilayer generalization of the IQH 
state. The topological properties of such systems are 
characterized by three spectral Chern numbers42i In our 
case, the nonzero topological invariant is given by the 
winding number of vector field m{'k±;kz) which maps 
the in-plane hexagonal Brillouin zone (a 2-torus) to a 
unit 2-sphere: 



If K 

An J y ^ ^ ' \k\ 



(15) 



By treating as a parameter, Hamiltonian ([5]) has ex- 
actly the same form as that of the 2D Kitaev model. 
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The first Chern number (|15l) of tlie corresponding '2D' 
Hamiltonian Hk^{'k±) is given hy ly = sgnA/^^)^ where 
Afe^ = 6^/3{k + Kc sinkz) is the effective gap parameter 
at corner K of the hexagonal Brillouin zone. Therefore, 
as long as |k| > Kc, hence the system remains gapped, 
the winding number v — sgnK = ±1 is independent of 
kz, and the whole Brillouin zone is characterized by the 
same chirality. 



IV. STAGGERED MAGNETIC FIELD 

The gapped phase discussed in the previous section is 
topologically trivial due to the absence of TRS. As dis- 
cussed in Ref. :18:, TRS is a prerequisite for the existence 
of 3D topological insulators. In this section, we show 
that one can introduce a gap to the fermion spectrum, 
while at the same time preserving an effective TRS, by 
subjecting the system to a staggered magnetic field, i.e. 
rjj — (— Because the sign of the field alternates be- 
tween successive honeycomb layers, a discrete symmetry 
emerges in our model system as Hamiltonian ([5]) is invari- 
ant under time-reversal T followed by a lattice translation 
along z axis. This additional symmetry manifests itself as 
a TRS in the auxiliary Majorana hopping problem. The 
mechanism proposed here is similar to Haldane's model 
of realizing 2D quantum Hall effect without a net mag- 
netic fiux through the unit cells4i 

Due to the staggered field, the sign of second nearest- 
neighbor hopping Uji^ also alternates between successive 
honeycomb planes. This results in a staggered gap func- 
tion ±A(kj^), and a doubled unit cell along z direction. 
We denote the fermion annihilation operators on the even 
and odd layers as a^.s and 6r,s, respectively, where sub- 
script s = 1 , 2 refers to the two inequivalent sites within a 
honeycomb plane. The Fourier-transformed Hamiltonian 
then reads H = ^ Ek *k-f^(k)*k, with 



Hik) 



/ A(ki) g{kz) -if{k^) 

g{kz) -A(ki) -«/(ki) 

ifik^r -A(ki) -gikz) 

V ^/(ki)* -gikz) A(ki) 



(16) 



= A(k_L)rV^+.9(fc,)T^CT^+Im/(k_L)r^+Rc/(k_L)T^, 

and = {au,i,b^,i,o,k,2,b^,2) ■ The two sets of Pauli 
matrices and a^^ now act on the sublattice (1,2) and 
even-odd (a, b) indices, respectively. 

In addition to symmetry relation (jlOp shared by Hamil- 
tonians describing free Majorana fermions, the hermitian 
matrix (ITBl) also satisfies 



T^iia^) H^{k) i-iay) = H{-k), 



(17) 



stemming from the generalized TRS. The extra factor 
can be gauged away by a 7r/2 rotation about the axis. 
Eq. (1171) defines the symmetry property of Dili Hamilto- 
nians in Altland-Zirnbauer's classification.'^*' As discussed 
in Ref. [l^ a topological invariant can be defined for 



Hamiltonians in this symmetry class based on the block 
off-diagonal representation of the Hamiltonian, or more 
precisely, of the spectral projection operator. In fact, the 
same definition can be applied to classes of Hamiltoni- 
ans which possess some form of chiral symmetry arising 
from either the sublattice symmetry or a combination of 
particle-hole and time-reversal symmetries.^* 

To compute the topological winding number of our 
model, we first bring Hamiltonian (jl6p into a block off- 
diagonal form through a series of unitary transforma- 
tions. First, noting that the layered honeycomb lattice 
is bipartite in which nearest neighbors of one sublattice 
belong to the other one, we regroup fermions of the same 
sublattice into a block, e.g. ak,i and 6k, 2- Mathemati- 
cally this is achieved by interchanging the 2nd and 4th 
entries of ^'k, the transformed Hamiltonian becomes 

H{k) ^ a^Im/(ki) + a^Re/(ki) + a'-g{kz) + /3A(ki), (18) 

where a'' and /3 given by 

= 7^ = r^cr'', ^ = 70 = r^ (19) 

are the standard Dirac matrices. The second part of the 
unitary transformation is a 7r/2 rotation about the new 
axis, which transforms — !■ r^, hence 



where the upper right block is 



D{k) 



g(fc,)-zA(ki) -if{k^) 
z/(ki)* -g{kz) -iA{k^) 



(20) 



(21) 



Noting that /(-k_L) = /(k_L)*, and A(k_L), g{kz) are odd 
functions of k, matrix (|21l) satisfies a symmetry (k) — 
—D{—k). The block off-diagonal form of the Hamiltonian 
implies that e(k)2 = D''{k)D{k), which gives rise to a 
fermion spectrum 



e±(k)=±VA(ki)2+g(fc, 



|/(k±)| 



(22) 



Due to the effective TRS (IT7|) . the spectrum is double de- 
generate at each wavevector k, except at possible Fermi 
points. It is worth noting that, contrary to the uniform 
field case, the staggered field immediately opens an en- 
ergy gap to the spectrum in the weak-pairing phase of 
the model. 

The topological properties of the quantum ground 
state (occupied Bloch states) in the gapped phase is cap- 
tured by the following Q matrix:^* 



0(k) = 





zt(k) 







'?(k) = 



^(k) 
£+(k)- 



(23) 



In fact Q(k) represents a 'simplified' Hamiltonian ob- 
tained by assigning an energy e = —1 to all occupied 
states and e = -1-1 to all empty bands of Hamiltonian 
ij'(k) ]^*'^* As long as the system is in the same gapped 
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(a) 



(b) 



K 



V=0^ 



, v=+l 



V=0 



V=0 



i 




v=+l 






v=o 


v=-l 





(c) 



(d) 



1 (K>^>Q-eH^-G-G-G-£^-- 



V 



OOOOOOOOO-O 



1 O-e-G-G-G-O-G-e-G-G-- 



—Q-Q-e-o-o-Q-e-e-3-^i 1 



( FG-G-G -G -Q- G- G-O-Q - - 



V 



0.02 0.04 
K 



FIG. 4: (Color online) (a) Phase diagram in the plane of 
Jx + Jy + Jz = const. The shaded triangle corresponds to 
the region where the in-plane couplings satisfy the triangle 
inequalities. In the absence of staggered field k = 0, the 
shaded area represents a critical phase with zero energy gap. 
We define a parameter S — 2Jx ~ Jz along the Jx ~ Jy line 
specifying the distance from the phase boundary, (b) Phase 
diagram as a function of distance 5 and field strength k. Pan- 
els (c) and (d) Numerical evaluation of the winding number 
as a function of S and k, respectively. The field strength 
K — ±0.02Jz in the calculation of panel (c), whereas we set 
Jx ~ Jy ~ Jz = J and K is measured in units of J in (d). 



phase, one can continuously deform the model such that 
i/(k) gradually transforms to the simplified form 
Not surprisingly, the Q matrix is related to the spectral 
operator via Q(k) = 1 - 2P{k)}^ 

It is easy to see that the block matrix q satisfies 



gt(k)<7(k) = l, (k) = -q{-k), 



(24) 



as expected for a Dili class Hamiltonian. The topological 
invariant characterizing the ground state is defined as the 
integer winding number of mapping q : — >■ U{2)^ii^ 



247r2 



l^''nT[{q-^d^q){q-^d.q){q-^dpq)\, (25) 



where the integral is over the three-dimensional Brillouin 
zone, which is essentially a 3-torus . For Dili class 
Hamiltonians, the winding number v can take on an ar- 
bitrary integer, each labels a unique topological class of 
the quantum ground state. 

The value of the topological winding number v can 
only change in the presence of a quantum phase transi- 
tion, which is usually accompanied by a vanishing bulk 
gap. In our case, noting that 17(0) = 0, a prerequisite for 
spectrum ([22)1 to be gapless is that the in-plane couplings 
satisfy the triangle inequalities | J^j < \Jx\ + | and so 
on [see Fig.lU^a)], such that solutions exist for /(k_L) = 
(the weak-paring regime of the model). In the absence 



of magnetic field k = /i = 0, the triangle inequalities 
thus define a critical phase with gapless fermionic excita- 
tions. For convenience, we use 5 to denote the 'distance' 
from the boundary of the critical phase. For example 
5 = 2Jx — Jz along the J^ = Jy line shown in Fig. HJ^a). 
We numerically compute the winding number v using 
definition ((25|) for various gapped phases of the model 
Hamiltonian; the resulting phase diagram as a function 
of 5 and k is summarized in Fig. lUJb). 

At the phase boundary defined by k = and (5 = 0, 
the system undergoes a quantum phase transition of a 
topological nature. Fig. Sl^c) shows numerical evaluation 
of the winding number as a function of ^ = 2 J^, — Jz in 
the presence of a staggered field such that |k| — 0.027^. 
Depending on the sign of k, the winding number jumps 
from z^ = Otoz^ = ±l when S crosses the phase boundary 
from the topologically trivial phase (corresponding to the 
A phase in 2D Kitaev model). On the other hand, for 
systems inside the critical phase (S < 0), the winding 
number jumps from zy = — ltoiy = +lasK changes sign. 

The nontrivial winding number of the quantum ground 
state can also be understood from the topological prop- 
erties of the singular (Dirac) points in the fermion spec- 
trum. To simplify the discussion, we focus on the sym- 



metric case Jx 
(1221) in the n - 



Jy Jz 



J, where zeros of spectrum 



limit are at 
kt = (f,0,0), k* = (-f,0,0). 



(26) 



In the vicinity of these two points, the fermions obey a 
Dirac-like Hamiltonian [c.f. Eq. (fTS)) ] 

H{k* + p) = a^p^ + m/3, (27) 

where the scaled momentum and mass term are 

p_L = V3J{Py,TPx), Pz = 4:Jj_Pz, 

m = A(k*) ±6V3k. 



(28) 



The plus and minus signs in the expression of m corre- 
spond to k^ and kg, respectively. The spectrum of (P7)) 
given by e± = ■iz^Jp^ + m? is two-fold degenerate at each 
p. As discussed above, a 7r/2 rotation about brings 
the Dirac mass into a chiral mass term, i.e. /3 — 17^/?. 
The transformed Hamiltonian H = a^p^ — if]j^m is in 
a block off-diagonal form with D = p ■ cr — im. The q 
matrix of the spectral projector is then given by 



g(k* + p) = 



p • (T — tm 



VP 



(29) 



Eqs. (P7)) and ([^^ can be viewed as a continuum approx- 
imation to the low-energy physics of the model system. A 
direct evaluation using Eq. (j25l) yields a winding number 



± 



1 m 

2 M' 



(30) 



with ± sign refering to k^ and kj points, respectively. 
Note that when evaluating v using the continuum de- 
scription, we have extended the domain of integration to 
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a 3-sphere in Eq. (|25|) . The appearance of a half- integer 
V is an artifact of the continuum description; the winding 
number is modified once contributions from high-energy 
Bloch states (away from the Dirac point in the Brihouin 
zone) are properly included. Since the mass term has 
opposite sign at the two Dirac points, r7i(k|) — —■m{V^2)^ 
we obtain v\ — ~ ^sgnK. Interestingly, the sum of 
these two winding numbers v — v\ + V2 reproduces the 
topological invariant of the lattice model. 

As recently pointed out in Ref. stable Fermi lines 
generally appears in 3D topological superconductors de- 
scribed by Hamiltonians belonging to classes CI or DHL 
An extended phase diagram of a lattice CI model^" in- 
deed shows regions of gapless phase with topologically 
stable Fermi lines-^- Here we show that such Fermi lines 
are also possible in our Dili Hamiltonian p^ . 

To this end, we consider perturbations which break 
the sublattice symmetry by introducing different second- 
neighbor hoppings K on the two inequivalent sites of the 
lattice: ki^2 = k± <^k. In the vicinity of the Dirac points 
(j26p . the asymmetric coupling results in an additional 
term of the off-diagonal block matrix 

Z? = p • cr — im — ifia^ , /i — 6V^6k. (31) 

The corresponding spectrum ^ fj'^ + (^fi— ^rn? -f p'^)^ 
has a Fermi line specified by 

Vi_ = v/Ai'-m2, ^ 0, (32) 

when fi > m. Finally, we point out that introducing 
a similar asymmetric hopping Sk in the case of uniform 
field results in a Fermi surface in the weak-pairing regime. 



V. GAPLESS SURFACE MAJORANA 
FERMIONS 

Analogous to the gapless chiral edge modes in the non- 
Abelian phase of 2D Kitaev modelj^ the nontrivial topol- 
ogy of the — ±1 quantum ground state in our 3D 
model manifests itself through the appearance of gap- 
less Majorana fermions at the sample surface. To study 
the properties of these surface modes, we solve the Ma- 
jorana hopping problem (j4]) and ([71) in a finite geometry 
with periodic boundary conditions along x and z direc- 
tions and open boundary condition along y direction. We 
assume that the x-axis is parallel to the zigzag direction 
of the honeycomb lattice. As Fig. [5] shows, in addition to 
gapped states in the bulk, the spectrum of a finite system 
contains additional surface modes crossing the bulk gap. 

The properties of these surface states can be under- 
stood analytically using the edge modes of 2D Kitaev 
model, whose existence has been demonstrated in Ref. [26l. 
These boundary modes are confined to the edge of the 
honeycomb lattice and possess a definite chirality de- 
pending on the sign of magnetic field h. Without loss 
of generality, we assume k ~ X^h > 0, hence e±{kx) > 



(a) k^ = (b) k„ = 7c 




It 27[ -7t/2 7c/2 



FIG. 5: (Color online) Spectrum of a finite system with L = 
10 layers along y direction. We consider the symmetric case 
with Jx ~ Jy = Jz = J- Coupling constant along vertical 
links is set to Jy — 0.35J, and the strength of the uniform 
field is K = 0.05J. Panels (a) and (b) show dispersion along 
kx and fez directions, respectively, of the 2D surface Brillouin 
zone. 

for positive in the ID Brillouin zone. The Hamiltonian 
describing the edge states is given by 

■Hedge = ^ X! Xi-kx) xikx), (33) 

where xi^x) and x{~kx) = xH^x) for kx > are anni- 
hilation and creation operators, respectively, of the Ma- 
jorana edge modes. The spectrum of the edge states has 
the form e±{kx) ~ 12k sin fc^; in the vicinity of the ID 
Fermi point k* — tt, and gradually merges with the bulk 
spectrum as kx moves away from k*^ 

It is interesting to note that the Jy — limit of the 3D 
model HI) can be viewed as a collection of decoupled hon- 
eycomb layers. Due to the staggered magnetic field, the 
edge states on even and odd numbered layers have oppo- 
site chirality v = ztsguK; the corresponding quasiparticle 
operators are denoted as x+(^a;,2n) and X-{kx,'2n + 1), 
respectively. The ensemble of the edge modes is described 
by Hamiltonian: 

'H± = ^^^ej_{k^)[x+{-k^,2n)x+{kx,2n) 

n 

^X-{-kx,2n + l)x-(k^,2n+l)]. (34) 

For odd-numbered layers, the quasiparticle annihilation 
operator is given by X-{kx, 2n + 1) with negative k^- A 
nonzero Jy introduces coupling between adjacent layers: 

^11 = ^^11 [x-i-kx,2n + l)x+{kx,2n) 

n 

+X+{~kx,2n)x-{kx,2n~\)\. (35) 

After Fourier transformation with respect to z coordi- 
nate, we obtain the Hamiltonian for surface modes 

Hsurf = ^ Y V'L,fe,[e-L(fc.)tT^ + e||(A:.)a"]Vfc.,fc.,(36) 

where ipk^M. = {x+{kx,kz),X-{kx,kz))'^ , and the off- 
diagonal coupling £11(^2) = 2J|| sinfc^. The spectrum can 
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1 /2 

be easily obtained esurt — (*^| + ^1) • Consistent with 
the numerical calculation shown in Fig. [SJ the surface 
spectrum has a conic singularity at the 2D Fermi point 
(fcj,fc*) = (7r,0). 

In the continuum approximation, the low-energy sur- 
face modes close to the 2D Fermi point obeys a gapless 
Dirac-like Hamiltonian -ffoirac = —i(u±CF^dx + u\\(T^dz)- 
The existence of gapless surface Majorana modes is a 
direct consequence of the nontrivial ground-state topol- 
ogy when the 3D bulk sample is terminated by a 2D 
boundary.>Si^ Due to its topological nature, the surface 
Dirac cone is stable against perturbations respecting the 
effective TRS. 



VI. DISCUSSION 

To summarize, we have constructed an interacting 
bosonic model on a three-dimensional layered honeycomb 
lattice which displays topologically nontrivial ground 
states. Our approach is based on the F-matrix general- 
ization of Kitaev's original spin-1/2 modeli^"— We show 
that the ground state in the weak-pairing phase remains 
gapless in the presence of a small magnetic field. The 
fermion spectrum of this critical phase is characterized 
by four topologically protected Fermi points. Upon in- 
creasing the field strength, Fermi points with opposite 
winding numbers move toward each other and eventually 
annihilate at a critical field, signaling the transition into 
a topologically trivial phase with a gapped spectrum. 

On the other hand, with the sign of magnetic field 
staggered along successive layers, the bulk spectrum of 
the weak-pairing phase immediately acquires a gap. An 



effective TRS is restored thanks to the staggering of the 
magnetic field. The corresponding quantum ground state 
is characterized by an integer winding invariant and pos- 
sesses nontrivial topological properties, a manifestation 
of which is the appearance of protected gapless surface 
Majorana modes when the sample is terminated by a 
two-dimensional surface. 

It is worth noting that our model provides an al- 
ternative example in which a 3D topological insulator 
emerges from an interacting bosonic Hamiltonian. More 
importantly, the TRS essential to the existence of 3D 
topological insulators is realized in our model through 
a generic physical mechanism, instead of relying on 
special symmetries of the F-matrix representation in 
the recently proposed diamond-lattice modeli^i^ We 
also remark that despite the experimental difficulty of 
realizing the Kitaev model and its variants, the study of 
these exactly solvable models has enriched our under- 
standing of the physics of topological phases. Finally, we 
would like to point out that it remains unclear whether 
an emergent topological insulator can be realized in a 
spin-1/2 Kitaev model on a 3D trivalent lattice. 
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